library(tidybulk)
library(tidySummarizedExperiment)
library(ggpubr)
library(GEOquery)
library(readxl)
library(tidyverse)

# bradley2018 --------
bradley <- read_csv('mission/HMCES_HIV/bradley2018_GSE115449_EffCountFile.csv.gz')

bradley <- bradley |>
  pivot_longer(2:last_col())

bra_meta <- read_delim('mission/HMCES_HIV/bradley2018.meta.txt',
                       col_names = c('acc','sample')) |>
  separate(sample, into = c('group','name'))

bradley <- bra_meta |>
  select(-acc) |>
  right_join(bradley)

bradley <- bradley$Transcript |>
  clusterProfiler::bitr(fromType = 'REFSEQ',
                        toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  rename(Transcript = REFSEQ) |>
  left_join(bradley) |>
  as_tibble()

tdb_bra <- bradley |>
  tidybulk(.sample = name, .transcript = SYMBOL,
           .abundance = value)

tdb_bra <- tdb_bra |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

tdb_bra |>
  ggplot(aes(value_scaled, group = name, color = group)) +
  geom_density() +
  scale_x_log10()

tdb_bra |>
  reduce_dimensions(method = 'PCA') |>
  pivot_sample() |>
  ggplot(aes(PC1,PC2, color = group)) +
  geom_point()

test_bra <- tdb_bra |>
  test_differential_abundance(~0 + group,
                              contrasts = c('groupBNab-groupControl'))

bra_sig <- test_bra |>
  keep_abundant() |>
  pivot_transcript()

bra_sig |>
  filter(SYMBOL == 'HMCES') |>
  select(`FDR___groupBNab-groupControl`)

tdb_bra |>
  filter(SYMBOL == 'HMCES') |>
  count(group)

tdb_bra |>
  filter(SYMBOL == 'HMCES') |>
  mutate(group = fct_relevel(group, 'Control')) |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 3000, yend = 3000) +
  annotate(geom = 'text',x = 1.5, y=3150, size =5,
           label = 'logFC=0.135, P.adj=0.0734') +
  scale_color_discrete(label = c('Ctrl (n=46)', 'BNab (n=46)')) +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in PBMC',
       subtitle = 'Bradley 2018 GSE115449')

# austin2019 -------
austin <- read_delim('mission/HMCES_HIV/Austin2019_GSE119234_LN_BcellSubset_raw_counts.txt.gz')

tdb_aus <- austin |>
  pivot_longer(2:last_col()) |>
  rename(symbol = ...1) |>
  mutate(sample = name) |>
  separate_wider_delim(name, delim = ':',
                       names = c('individual','cell_type')) |>
  tidybulk(.sample = sample,
           .transcript = symbol,
           .abundance = value)

tdb_aus <- tdb_aus |>
  identify_abundant(factor_of_interest = individual) |>
  scale_abundance()

tdb_aus |>
  ggplot(aes(value_scaled, group = sample, color = individual)) +
  geom_density() +
  scale_x_log10()

tdb_aus <- tdb_aus |>
  mutate(cell_type = str_remove(cell_type, '_HIV'))

tdb_aus <- tdb_aus |>
  mutate(infect = str_extract(individual, 'HD|HIV'))

tdb_aus |>
  reduce_dimensions(method = 'PCA') |>
  pivot_sample() |>
  ggplot(aes(PC1,PC2, color = cell_type, shape = infect)) +
  geom_point()

tdb_aus |>
  filter(symbol == 'C3orf37') |>
  ggplot(aes(individual, cell_type, fill = log1p(value_scaled))) +
  geom_raster() +
  scale_fill_gradient(high = 'red', low = 'blue') +
  theme_pubr(legend = 'right') +
  labs(fill = 'log-Normalized expression',
       title = 'HMCES expression in B cell subsets',
       subtitle = 'Austin 2019 GSE119234')

tdb_aus |>
  filter(symbol == 'C3orf37') |>
  as_tibble() |>
  tidyHeatmap::heatmap(individual, cell_type, value_scaled,transform = log1p,
                       scale = 'column',
                       palette_value = c('blue','grey','red'))

## compare HIV GC vs Ctrl GC --------
aus_hivhd <- tdb_aus |>
  filter(cell_type == 'GC') |>
  test_differential_abundance(~0+infect,
                              contrasts = 'infectHIV-infectHD')

aus_hivhd |>
  pivot_transcript() |>
  filter(symbol == 'C3orf37') |>
  select(`FDR___infectHIV-infectHD`)

aus_hivhd |>
  filter(symbol == 'C3orf37') |>
  ggplot(aes(infect, value_scaled, color = infect)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .01) +
  annotate(geom = 'segment',x = 1, xend = 2, y=6100, yend =6100) +
  annotate(geom = 'text',x = 1.5, y=6300, size =5,
           label = 'logFC=-0.213, P.adj=0.65') +
  theme_pubr() +
  labs(y = 'Normalized expression', color = 'group',
       title = 'HMCES expression in GC B cell',
       subtitle = 'Austin 2019 GSE119234') +
  scale_color_manual(values = c('blue','red'),
                     label = c('HD (n=5)','HIV (n=4)'))

aus_hivhd |>
  filter(symbol == 'C3orf37' & infect == 'HIV')

## BNab-hi vs BNab-lo -----------
aus_bnab <- tdb_aus |>
  mutate(bnab = case_when(individual %in% c('HIV1','HIV4') ~ 'high',
                          infect == 'HIV' ~ 'low',
                          .default = NA)) |>
  filter(cell_type == 'GC' & infect == 'HIV') |>
  test_differential_abundance(~0+bnab,
                              contrasts = 'bnabhigh-bnablow')

aus_bnab |>
  pivot_transcript() |>
  filter(symbol == 'C3orf37') |>
  select(`logFC___bnabhigh-bnablow`,`FDR___bnabhigh-bnablow`)

aus_bnab |>
  filter(symbol == 'C3orf37') |>
  mutate(bnab = fct_relevel(bnab, 'low')) |>
  ggplot(aes(bnab, value_scaled, color = bnab)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .01) +
  annotate(geom = 'segment',x = 1, xend = 2, y=6150, yend =6150) +
  annotate(geom = 'text',x = 1.5, y=6350, size =5,
           label = 'logFC=0.505, P.adj=0.992') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in GC B cell',
       subtitle = 'Austin 2019 GSE119234') +
  scale_color_manual(values = c('red','blue'))

# martin2020 ---------
meta_mart <- getGEO('GSE141498')

meta_mart <-
meta_mart$GSE141498_series_matrix.txt.gz@phenoData@data |>
  select(title, `cohort:ch1`) |>
  as_tibble()

meta_mart <- meta_mart |>
  mutate(name = str_replace(title, ' for Patient ', '_'),
         id = str_extract(name, 'EMG.+'),
         cell = str_extract(name, '.+(?=_EMG)'),
         cohort = `cohort:ch1`,.keep = 'none')

martin <- read_csv('mission/HMCES_HIV/martin2020_GSE141498_RawCounts.csv.gz')

m <- martin |>
  separate_wider_delim(ID, names = c('ENSEMBL','SYMBOL'),delim = '|') |>
  select(-ENSEMBL) |>
  pivot_longer(-1) |>
  left_join(meta_mart)

tdb_mar <- m |>
  tidybulk(.sample = name,
           .transcript = SYMBOL,
           .abundance = value)

tdb_mar <- tdb_mar |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = cohort) |>
  scale_abundance()

tdb_mar |>
  ggplot(aes(value_scaled, group = id, color = cohort)) +
  geom_density() +
  scale_x_log10()

tdb_mar |>
  reduce_dimensions(method = 'PCA') |>
  pivot_sample() |>
  ggplot(aes(PC1,PC2, color = cohort, shape = cell)) +
  geom_point()

tdb_mar |>
  filter(SYMBOL == 'HMCES') |>
  ggplot(aes(cell, id, fill = log1p(value_scaled))) +
  geom_raster() +
  scale_fill_gradient(high = 'red', low = 'blue') +
  theme_pubr(legend = 'right') +
  labs(fill = 'Normalized expression')

## compare Bnab vs no_Bnab --------
mar_neu <- tdb_mar |>
  filter(cell == 'Bcells') |>
  test_differential_abundance(~0+cohort,
                              contrasts = 'cohortNeut-cohortNonNeut')

mar_neu |>
  pivot_transcript() |>
  filter(SYMBOL == 'HMCES') |>
  select(`logFC___cohortNeut-cohortNonNeut`,`FDR___cohortNeut-cohortNonNeut`)


tdb_mar |>
  filter(cell == 'monocytes') |>
  test_differential_abundance(~0+cohort,
                              contrasts = 'cohortNeut-cohortNonNeut') |>
  pivot_transcript() |>
  filter(SYMBOL == 'HMCES') |>
  select(`logFC___cohortNeut-cohortNonNeut`,`FDR___cohortNeut-cohortNonNeut`)

mar_neu |>
  filter(SYMBOL == 'HMCES') |>
  mutate(cohort = fct_relevel(cohort, 'NonNeut')) |>
  ggplot(aes(cohort, value_scaled, color = cohort)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 750, yend = 750) +
  annotate(geom = 'text',x = 1.5, y=800, size =5,
           label = 'logFC=-0.498, P.adj=0.721') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in B cell',
       subtitle = 'Martin 2020 GSE141498') +
  scale_color_manual(values = c('red','blue'),
                     label = c('No-bnAb (n=35)','bnAb (n=14)'))

tdb_mar |>
  filter(SYMBOL == 'HMCES' & cell != 'Bcells') |>
  mutate(cohort = fct_relevel(cohort, 'NonNeut')) |>
  ggplot(aes(cohort, value_scaled, color = cohort)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 750, yend = 750) +
  annotate(geom = 'text',x = 1.5, y=800, size =5,
           label = 'NS') +
  theme_pubr() +
  facet_wrap(~cell) +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in other immune cells',
       subtitle = 'Martin 2020 GSE141498') +
  scale_color_manual(values = c('red','blue'),
                     label = c('No-bnAb (n=35)','bnAb (n=14)'))

# luo2022 ----------
luo22 <- read_excel('mission/HMCES_HIV/luo2022_GSE157198_readcount_GEO.xlsx')

luo_meta <- read_delim('mission/HMCES_HIV/luo2022_meta.txt',
                       col_names = c('acc','sample'))

luo_meta <- luo_meta |>
  mutate(name = str_extract(sample, 'S\\d+'),
         group = str_extract(sample, 'EC|Healthy|ART-na|ART') |>
           make.names(), .keep = 'none')

tdb_luo <- luo22$geneID |>
  clusterProfiler::bitr(fromType = 'ENSEMBL',
                        toType = 'SYMBOL',
                        OrgDb = 'org.Hs.eg.db') |>
  rename(geneID = ENSEMBL) |>
  left_join(luo22) |>
  select(-geneID) |>
  as_tibble() |>
  pivot_longer(-1) |>
  left_join(luo_meta)

# identify abundance transcripts is highly recommended before scaling
tdb_luo <- tdb_luo |>
  tidybulk(.sample = name,
           .transcript = SYMBOL,
           .abundance = value) |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

tdb_luo |>
  reduce_dimensions(method = 'pca') |>
  pivot_sample() |>
  ggplot(aes(PC1, PC2, color = group)) +
  geom_point()

tdb_luo |>
  mutate(group = make.names(group)) |>
  test_differential_abundance(~0+group,
                              contrasts = 'groupEC-groupART.na') |>
  filter(SYMBOL == 'HMCES') |>
  as.data.frame()

tdb_luo |>
  filter(SYMBOL == 'HMCES') |>
  ggplot(aes(group, value_scaled, color = group)) +
  stat_mean(geom = 'col', fill = 'white') +
  geom_jitter(height = 0, width = .1) +
  #annotate(geom = 'segment', x = 1, xend = 2, y = 750, yend = 750) +
  #annotate(geom = 'text',x = 1.5, y=800, size =5,
  #         label = 'NS') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in monocytes',
       subtitle = 'Luo 2022 GSE157198')

# parker2023 ----------
parker <- data.table::fread('mission/HMCES_HIV/parker2023_GSE199911_processeddata_genecounts.csv.gz')

parker[1:5,1:8]

parker <- parker |>
  mutate(group = ifelse(studygroup == 1, 'HIV', 'HC')) |>
  select(-c(V1,studygroup:age)) |>
  pivot_longer(where(is.numeric)) |>
  mutate(ENSEMBL = str_remove(name, '\\..+')) |>
  select(-name)

parker <- parker |>
  ensembl_to_symbol(.ensembl = ENSEMBL) |>
  filter(!is.na(transcript))

tdb_par <- parker |>
  as_tibble() |>
  tidybulk(.sample = sample,
           .transcript = SYMBOL,
           .abundance = value) |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group) |>
  scale_abundance()

tdb_par |>
  reduce_dimensions(method = 'pca') |>
  pivot_sample() |>
  ggplot(aes(PC1, PC2, color = group)) +
  geom_point()

tdb_par |>
  test_differential_abundance(~0+group,
                              contrasts = 'groupHIV-groupHC') |>
  pivot_transcript() |>
  filter(SYMBOL == 'HMCES') |>
  as.data.frame()

tdb_par |>
  filter(SYMBOL == 'HMCES') |>
  ggplot(aes(group, value_scaled, color = group)) +
  geom_boxplot() +
  geom_jitter(height = 0, width = .1) +
  annotate(geom = 'segment', x = 1, xend = 2, y = 3100, yend = 3100) +
  annotate(geom = 'text',x = 1.5, y=3200, size =5,
           label = 'NS') +
  theme_pubr() +
  labs(y = 'Normalized expression',
       title = 'HMCES expression in PBMC',
       subtitle = 'Parker 2023 GSE199911')
